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摘要 


本 文 基于 阻尼 块 反 帘 法 与 子 空间 投影 算法 设计 了 一 种 求解 特征 值 问 题 的 广义 
共 轿 梯度 算法 , 同时 也 实现 了 相应 的 计算 软件 包 ， 然 后 对 算法 和 计算 过 程 进行 了 
一 系列 的 优化 来 提高 算法 的 稳定 性 、 计 算 效 率 和 并 行 可 扩展 性 , 使 得 本 文 的 算法 适 
合 在 并 行 计 算 环 境 下 求解 大 规模 稀 朴 矩阵 的 特征 值 . 所 形成 的 软件 包 是 基于 Matrix- 
Free 和 Vector-Free 设 计 的 , 可 以 应 用 于 任意 的 矩阵 向 量 结构 .针对 几 种 典型 算 阵 的 测 
试 结果 表明 本 文 的 算法 和 软件 包 不 但 具有 良好 的 数值 稳定 性 , 同时 相 比 于 SLEPc 软 
件 包 中 的 LOBPCG 以 及 Jacobi-Davidson 解 法 器 有 2-6 倍 的 效率 提升 ， 软 件 包 的 网 址 : 
https://github.com/pase2017/GCGE-1.0. 
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A Generalized Conjugate Gradient Method for 
Eigenvalue Problems 


Abstract 


A generalized conjugate gradient method is proposed to solve eigenvalue problems. 
This method is designed by combining the dumping block inverse power scheme, subspace 
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projection method. Furthermore, based on the properties of the proposed method, a series 
of optimization techniques is developed to improve the stability, computing efficiency and 
scalability. We also introduce a computing package GCGE (Generalized Conjugate Gradi- 
ent Eigensolver) which is developed based on the proposed method here. Some numerical 
examples are provided to validate the stability, computing efficiency and scalability of the 
method in this paper. The corresponding computing package can be downloaded from 
the web site: https://github.com/pase2017/GCGE-1.0. 


Keyword: eigenvalue problem; dumping block inverse power; generalized conjugate 
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1 算法 介绍 


本 文 介 绍 一 种 基于 阻尼 块 反 寺 法 的 并 行 特征 值 求解 算法 及 其 相应 的 高 效 实现 方法 . 
阻尼 块 反 圭 法 是 反 堪 法、 阻尼 思想 和 子 空 间 投影 方法 的 组 合 , 同时 考虑 具体 的 并 行 实现 
效率 来 设计 出 相应 的 代数 特征 值 求解 算法 . 这 里 设计 的 算法 主要 是 用 来 求解 对 称 正定 矩 
阵 最 小 端 特征 值 及 其 相应 的 特征 向 量 . 当然 本 文 算法 设计 的 思想 也 可 以 用 来 构造 求解 非 
对 称 与 非 正定 矩阵 特征 值 问题 的 方法 . 在 具体 设计 的 时 候 主要 是 利用 阻尼 块 反 窜 法 结合 
并 行 计算 的 特点 来 设计 高 效 的 特征 值 求解 算法 . 

本 文 主要 考虑 求解 如 下 的 广义 代数 特征 值 问题 


Ax = ABr. 


为 了 简单 起 见 , 这 里 规定 矩阵 A4 和 B 是 对 称 正定 的 . 现 有 比较 流行 的 特征 值 解法 器 一 般 都 
基于 Krylov 子 空间 算法 [15], 比如 Arnoldi 算 法 、LOBPCG 算 法 和 Jacobi-Davidson 算 法 等 . 
基于 Krylov 子 空间 的 算法 (Arnoldi 算 法 和 Lanczos 算 法 |[15|) 的 特征 值 解法 器 , 需要 精确 求 
解 线 性 方程 组 , 否则 得 不 到 相应 的 上 Hessenberg 和 矩阵 .对 于 一 些 条 件数 很 大 的 矩阵 , 使 
用 迁 代 法 难以 收敛 , 而 使 用 直接 法 求解 , 会 使 得 求解 时 并 行 效率 较 差 . 基于 LOBPCG 或 
者 Jacobi-Davidson 等 算法 的 特征 值 解法 器 , 不 需要 精确 求解 线性 方程 组 , 因此 即使 对 于 
条 件数 很 大 的 矩阵 , 也 可 以 使 用 迭代 法 求解 , 并 行 效率 比较 好 . 但 目前 流行 的 特征 值 解 
法 器 中 , LOBPCG 稳 定性 较 差 , 这 是 由 其 产生 子 空间 的 方式 以 及 正 交 化 的 方式 所 导致 的 ， 
使 得 LOBPCG 在 求解 一 些 条 件数 较 大 或 者 精度 要 求 较 高 的 特征 值 问题 时 往往 求解 失败 . 
而 Jacobi-Davidson 算 法 需要 求解 不 定 且 几乎 奇异 的 线性 方程 组 , 需要 专门 设计 相应 的 线 
性 解法 器 来 进行 求解 . 在 使 用 Jacobi-Davidson 解 法 器 的 时 候 需 要 设置 的 参数 较 多 , 往往 
难以 选择 最 优 的 计算 参数 而 使 得 计算 效率 较 低 . 求解 大 规模 特征 值 问 题 的 并 行 算法 应 满 
足 如 下 条 件 : 一 方面 该 算法 不 需要 精确 求解 其 中 所 包含 的 线性 方程 组 , 另 一 方面 该 算法 
需要 具有 稳定 、 高 效 的 特点 才能 适合 求解 大 规模 矩阵 特征 值 问 题 的 要 求 . 同时 为 了 提高 
算法 的 并 行 可 扩展 性 , 尽量 把 算法 的 主要 部 分 转化 成 具有 良好 可 扩展 性 的 计算 方式 . 
基于 上 面 的 考虑 , 本 文 的 第 一 个 目的 是 设计 一 个 稳定 、 高 效 、 高 可 扩展 性 的 特征 值 
算法 . 为 了 适应 在 并 行 机 上 求解 大 规模 算 阵 的 特征 值 问 题 , 此 算法 不 需要 精确 求解 其 
中 所 包含 的 线性 方程 组 , 并 且 同 时 具有 和 稳定、 高效、 高 扩展 性 的 特点 . 我 们 基于 阻尼 块 
反 窜 法 和 子 空间 投影 算法 的 思想 设计 一 种 求解 特征 值 问 题 的 广义 共 罗 梯度 (Generalized 
Conjugate Gradient, 简称 GCG) 算 法 . 与 LOBPCG 算 法 相 比 较 , 主要 的 区 别 在 于 计算 向 
量 组 W 的 方式 , GCG 利 用 反 窜 法 结合 共 轮 梯度 (CG) 达 代 的 方式 , 而 LOBPCG 方 法 利用 的 
计算 残 差 向 量 组 的 方式 [5, 8, 11, 12, 13, 16]. 我 们 知道 随 着 特征 对 精度 的 增加 , 残 差 将 
变 得 越 来 越 小 从 而 使 得 LOBPCG 在 数值 上 不 稳定 5], 而 反 容 法 没有 这 种 缺点 . 当然 在 
没有 任何 机 器 误差 的 情况 下 这 两 种 方式 是 等 价 的 [5, 17]. 
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本 文 的 第 三 个 目的 是 基于 上 面 设计 的 特征 值 算法 及 其 性 质 设 计 一 个 具有 比较 好 的 
数值 稳定 性 (和 鲁 棒 性 ) 和 可 扩展 性 的 特征 值 并 行 解法 器 . 由 于 算法 主要 是 基于 矩阵 向 
量 的 整体 操作 来 运行 的 , 所 以 这 里 的 解法 器 可 以 设计 成 Matrix-Free 和 Vector-Free 的 形式 ， 
也 就 是 用 户 可 以 基于 自己 的 矩阵 和 向 量 格式 及 矩阵 和 向 量 操作 来 产生 具体 的 特征 值 解 
法 器 . 同时 为 了 提高 解法 器 的 效率 , 本 文 也 将 介绍 一 些 基 于 算法 性 质 的 优化 设计 . 

本 文 接 下 来 的 安排 为 : 本 节 主 要 介绍 GCG 算 法 , 第 2 节 介 绍 基 于 GCG 算 法 的 软件 包 ， 
第 3 节 介 绍 对 GCG 算 法 的 具体 实现 及 其 改进 , 第 4 节 将 进行 一 些 数值 实验 来 验证 本 文 提出 
算法 的 计算 效率 和 并 行 效率 . 最 后 一 节 给 出 对 本 文 的 总 结 . 


131 GCG 算 法 


GCG 算 法 是 一 种 子 空间 迭代 算法 , 它 用 阻尼 块 反 窜 法 达 代 的 方式 来 生成 一 个 三 元 问 
量 组 [X, P, W], 并 以 之 张 成 求解 特征 值 问 题 的 子 空间 . 三 元 同 量 组 |[X, P, WP XER 
次 迭代 中 的 近似 特征 向 量 , P 是 本 次 迭代 中 的 近似 特征 问 量 减 去 上 次 达 代 中 近似 特征 方 
向 的 分 量 , 丈 是 对 X 进 行 一 步 不 精确 反 里 法 和 迭代 得 到 的 向 量 . 由 于 向 量 组 中 的 W 是 通过 
— AN RG TP Be BI V AE RI, 对 于 对 称 正 定 问 题 , 我 们 常 采用 一 定 步 数 的 CG 和 迭代 来 
得 到 , Alc SSL AY SSA RE GUE (Generalized Conjugate Gradient Method). 
假设 要 求解 nev 个 特征 值 , 算法 1 给 出 了 具体 的 GCG 算 法 过 程 . 


Algorithm 1 GCG 算 法 

1. 初始 选取 XX, P= [], W =|], 其 中 XX 中 有 nev 个 向 量 . 计算 小 规模 特征 值 问题 
XTAXC = XTBXCOA 得 到 特征 值 A 和 特征 问 量 C. 更 新 X = XC. 

2. 对 线性 方程 组 W = A-1(BXA) 以 X 为 初始 值 执行 一 定 的 CG 迭代 步 得 到 向 
量 组 W. 

3. 定义 空间 V = [X, P, W], 对 其 进行 正 交 化 得 到 关于 和 矩阵 B 的 单位 正 交 基 组 V. 

4. 计算 Rayleigh-Ritz 问 题 : VTAVC = V7BVCA, 利用 得 到 的 特征 向 量 C 和 Ritz 
值 A 获 得 新 的 特征 向 量 逼 近 Xae。 = VC 和 相应 的 特征 值 通 近 . 

5. 检查 Xe。 的 收敛 性 , 如 果 收 敛 的 个 数 等 于 nev, 计算 结束 . 

6. 否则 , 计算 P = XQQVX, 更 新 X = Xrov, 回 到 步 2? 进 行 下 一 次 迭代 . 


注 1. 算法 1 步骤 6 中 的 Xuos\ 义 表示 在 Xes 中 去 掉 针 方向 上 的 分 量 , 后 文 有 对 具体 实现 方 
式 的 描述 . 

算法 1 与 LOBPCG 算 法 的 区 别 主要 在 于 生成 向 量 W 的 方式 和 对 向 量 组 V 的 正 交 化 
方式 . 算法 1 采用 块 反 徊 法 的 方式 来 产生 W, 同时 为 了 确保 算法 的 稳定 性 , 这 里 对 向 量 
组 V 进 行 完 全 的 正 交 化 . LOBPCG 方 法 中 的 向 量 组 全 是 当前 特征 通 近 的 残 差 向 量 , 并 
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HX, P, W 各 自 自 身 进行 正 交 化 , 并 不 对 整体 向 量 组 V 进 行 正 交 化 . 首先 当 近 似 特征 向 
量具 有 一 定 精 度 的 时 候 , 残 差 向 量 将 会 变 得 很 小 而 带 来 数值 不 稳定 性 . 另外 由 于 没有 
进行 整体 正 交 化 使 得 V 具 有 线性 相关 性 而 导致 在 计算 Rayleigh-Ritz 问 题 的 时 候 算法 终 
止 [5, 8]. 这 也 就 是 在 算法 1 中 我 们 要 对 向 量 组 V 进 行 正 交 化 的 原因 , 并 且 我 们 后 面 会 根据 
算法 的 性 质 来 提高 完全 正 交 化 的 计算 效率 . 

在 本 文 之 后 的 部 分 , 我 们 将 针对 算法 1 进行 一 系列 具体 实现 的 考虑 和 效率 的 优化 ， 
使 得 算法 在 保持 稳定 性 的 同时 具有 更 高 的 计算 效率 和 并 行 效率 , 以 及 可 以 进行 非 正定 
和 矩阵 特征 值 问题 的 计算 , 主要 包括 分 批 计算 技巧 、 块 计算 的 方式 、 块 正 交 化 、 快 速 计 
算 Rayleigh-Ritz 问 题 等 技术 , 其 目的 是 为 了 提高 计算 效率 和 算法 的 并 行 可 扩展 性 . 


1.2 ”特征 值 分 批 计算 


算法 1 中 求解 Rayleigh-Ritz 问 题 的 时 候 , 计算 特征 值 的 部 分 是 串 行 的 , 这 一 部 分 的 计 
算 时 间 很 难 通过 增加 进程 数 来 减少 . 当 计 算 的 特征 值 个 数 较 多 时 , 由 于 子 空间 特征 值 计 
算 时 间 与 特征 值 个 数 是 超 线性 关系 , 这 样 就 造成 总 的 计算 时 间 随 着 特征 值 个 数 的 增长 
是 超 线性 的 . 基于 这 样 的 考虑 , 算法 的 构造 应 该 尽量 减少 子 空间 特征 值 问题 的 维 数 . 2 
需要 计算 的 特征 值 个 数 较 多 时 , 我 们 将 特征 值 进行 分 批 计算 . EUCH VE SERA P I SCIL I 
AW fel. 算法 过 程 见 算 法 2. 
Algorithm 2 GCG 算 法 -特征 值 分 批 计算 
1. 初始 选取 [Xi Xe], Xo2[], P=[], W =|], 其 中 Xi 中 有 k 个 向 量 , Xo El 
含 已 经 收敛 的 特征 向 量 组 . 定义 向 量 组 X = [Xo, X4, Xo]. 计算 特征 值 问题 
XTAXC = XT BXCA. Wii, Xo] = XC. 
. 对 线性 方程 组 W = A! BX A) WX ABST — Ee I CGIK RAE S UI] 
量 组 W. 
. 定义 空间 V = [Xo, X1, X5, P,W], 并 计算 Rayleigh-Ritz 问 题 : VT AVC = 
VT BVOCN, 获得 新 的 特征 向 量 逼 近 [X?e”, X23"| 和 相应 的 特征 值 逼近 . 
4. 检查 收敛 性 : 更 新 [XX。, X1, Xo] = [X Xz**], Xo = [Xo, X], MAX WAI 
代 新 收敛 的 向 量 组 , Xi 为 未 收敛 向 量 的 前 k 个 , 其 余 为 XX2. iu X? 为 [X1, X] 
中 与 入 对 应 的 部 分 . 
5. 如 果 Xo 中 的 向 量 个 数 等 于 nev, 计算 结束 
. 否则 , 计算 P = X, — X24, iapx, Xa] = [Xi, Xo], 回 到 步 2 进行 下 一 次 迭代 . 


N 


CD 


c 


上 面 对 Rayleigh-Ritz 问 题 求 解 的 处 理 也 是 本 文 与 文章 回 中 的 区 别 之 一 , 这 里 只 需 
求解 Rayleigh-Ritz 问 题 的 目标 特征 对 , 而 不 需要 求解 全 部 的 特征 对 .通过 上 面 的 处 理 ， 
我 们 可 以 进一步 减少 求解 Rayleigh-Ritz 问 题 的 时 间 在 总 计算 时 间 中 所 占 的 比率 . 由 


于 Rayleigh-Ritz 问 题 是 一 个 稠密 矩阵 的 特征 值 问 题 , 并 行 计算 很 难 发 挥 出 计算 潜力 , 需 


要 特别 的 考虑 , 在 后 面 的 小 节 3.3 和 3.4 中 有 进一步 的 讨论 . 


1.3 ”市 位 移 的 QCG 算 法 


由 于 在 算法 2 的 第 2 步 中 计算 环 向 量 组 时 默认 使 用 CG 和 迭代 方法 , 而 CG 迭代 方法 只 适 


用 于 对 称 正定 线性 方程 组 的 求解 . 为 了 求解 对 称 不 定 和 矩阵 的 特征 值 问 题 , 我 们 对 算法 进 
行 了 位 移 处 理 , 使 得 算法 也 可 以 处 理 非 正定 的 问题 . 当然 对 于 更 广泛 的 情形 , 未 来 也 会 考 
FEMINRESIK (RATE AIGMRESIR (RAE. 


我 们 基于 算法 2 进行 修改 , 同样 考虑 对 特征 值 进行 分 批 计算 , 且 每 次 只 计算 K 个 已 向 量 


MAW HE. 算法 过 程 调 整 如 下 , 其 中 位 移 的 更 新 方式 表示 使 用 位 移 后 最 大 特征 值 为 最 
小 特征 值 的 100 倍 , 这 样 既 能 保证 特征 值 是 正 的 , 又 能 使 得 第 一 个 特征 值 不 会 太 小 . 


Algorithm 3 带 位 移 的 GCG 算 法 


1. 


初始 选取 [Xi, X5], Xo  [], P=[], W =[], 其 中 XX 中 有 k 个 向 量 , Xo 包含 已 
经 收敛 的 特征 向 量 组 . 定义 向 量 组 X = [Xo, Xi, Xo]. 令 位 移 参数 shift — 0. it 
算 特 征 值 问 题 X7AXC = XTBXCA. 更 新 [Xi1, X] = XC. 


. 对 线性 方程 组 W = (A + shift - B)-1(BXIA) 以 XX 为 初始 值 执 行 一 定 的 CG 适 


代步 得 到 新 的 向 量 组 W. 


. 定义 空间 V = [Xo, X4, Xo, P, W], 并 计算 Rayleigh-Ritz 问 题 : 


Vi(A+shift. B)VC= VT BVCA, 


获得 新 的 特征 向 量 逼 近 [X?s”, X3s 和 相应 的 特征 值 逼近 , 并 将 特征 值 减 去 
shift 以 得 到 原 问题 的 近似 特征 值 . 


. 检查 收敛 性 : HEBEL, X4, Xo] = [Xe X29], Xo = [Xo, X], deb X IVE 


代 新 收敛 的 向 量 组 , X NAUSEA AS, 其 余 为 X2. X AX, Xo] 
中 与 XX 对 应 的 部 分 . 更 新 shift = (nev — 100 - A1)/99. 


. 如 果 Xo 中 的 向 量 个 数 等 于 nev, 计算 结束 . 
. 否则 , 计算 已 = X, — X0, id, Xo] = [Xi Xo], 回 到 步 2 进行 下 一 次 迭代 . 


2 RETA 
基于 上 节 介 绍 的 GCG 算 法 , iH gii T REK IAE NP TE ISLI T" SCHERER 


度 解 法 器 (Generalized Conjugate Gradient Eigensolver, 简称 GCGE) 软 件 包 , 这 个 软件 包 


是 Matrix-Free 和 Vector-Free 的 . 在 该 软件 包 中 , 我 们 提供 了 基于 CSR. (B fr 47 HA Fi AR. 
阵 ), Hypre [10], PASE, PETSc [4], PHG [20], 以 及 SLEPc [9| 软 件 包 中 的 矩阵 结构 、 向 量 
结构 以 及 矩阵 向 量 操 作 的 特征 值 解法 器 的 调用 方式 . 如 果 用 户 使 用 的 是 其 中 一 种 矩阵 向 
量 结构 , 就 可 以 直接 使 用 软件 包 中 的 例子 程序 来 求解 特征 值 问 题 . 此 外 , 该 软件 包 支 持 用 
户 提 供 的 任意 矩阵 结构 、 疝 量 结构 以 及 秆 阵 问 量 操作 . 用 户 只 需 参 照 已 有 的 调用 方式 提 
供 自己 的 矩阵 结构 、 向 量 结构 、 和 矩阵 向 量 的 基本 操作 , 以 及 相应 的 接口 函数 , 就 可 以 使 
用 GCGE 软 件 包 在 自己 的 矩阵 癌 量 结构 下 求解 矩阵 特征 值 问题 . 具体 的 软件 包 C 语 言 代 
码 可 以 在 GitHub 上 找到 (网 址 :https://github.com/pase2017/GCGE-1.0), 欢迎 下 载 并 
测试 与 使 用 . 


2.1 程序 结构 


GCGE 软 件 包 主 要 分 为 gcge, app, example 三 个 部 分 . 其 中 , gcge 为 GCGE 的 主体 部 
分 , 用 来 具体 实现 GCG 算法 , 这 部 分 是 Matrix-Free 和 Vector-Free 的 , 不 涉及 任何 具体 的 
矩阵 结构 、 向 量 结构 与 矩阵 向 量 操作 ; app 提 供 有 具体 的 矩阵 向 量 结构 与 矩阵 向 量 操作 下 
对 gcge 的 调用 接口 , 目前 提供 了 CSR, Hypre, PASE, PETSc, PHG 以 及 SLEPc 等 软件 包 的 
调用 接口 ; example 中 给 出 了 相应 的 测试 和 调用 的 例子 . 

以 SLEPc 格 式 的 矩阵 向 量 结构 为 例 , 头 文 件 gcge.h 中 包含 了 JGCGE 算 法 求解 代数 特 
征 值 问题 的 所 有 调用 函数 以 及 参数 设置 函数 ; app 文 件 夹 中 的 头 文件 gcge_app_slepc.h 中 
提供 了 SLEPc 格 式 的 矩阵 向 量 结构 对 gcge 的 调用 接口 . 用 户 使 用 SLEPc 格 式 的 矩阵 向 量 
结构 计算 和 矩阵 特征 值 时 , 只 需要 包含 gcge.h 与 gcge_app_slepc.h 这 两 个 头 文件 就 可 以 调 
用 GCGE 软 件 包 中 提供 的 函数 进行 矩阵 特征 值 求解 的 操作 . 


2.1.1 ”内 核 结构 


头 文件 gcgeh 中 包含 了 GCG 算 法 求解 代数 特征 值 问题 的 所 有 调用 函数 以 及 参数 
设置 函数 .该 文件 中 包含 了 gcge_type.h，gcge_ops.h，gcge_para.h，gcge workspace.h, 
gcge rayleighritz.h, gcge orthogonalize.h, gcge xpw.h, gcge eigsolh, gcge solver.h^$ 
头 文件 , 他 们 共同 构成 了 实现 GCG 算 法 的 基础 . 

第 一 层 也 就 是 最 底层 的 头 文件 gcge_type.h 中 定义 了 一 些 基本 类 型 和 预定 义 , 这 个 头 
文件 会 被 gcge 中 的 所 有 其 他 头 文 件 所 包含 . 

第 二 层 含 有 两 个 头 文件 gcge_ops.h 和 gcge_Ppara.h.， 其 中 在 文件 gcge_ops.h 中 定义 
了 GCGE_OPS 结 构 体 , 该 结构 体 中 包含 了 需要 用 户 提 供 的 各 种 操作 . 这 其 中 包括 一 些 用 
户 必 须 提 供 的 基本 操作 , 以 及 用 户 可 以 选择 自己 提供 或 者 使 用 默认 操作 的 向 量 组 相关 的 
操作 ; 此 外 GCGE_OPS 中 还 提供 了 线性 求解 器 的 接口 , 用 户 可 以 选择 使 用 自己 提供 的 线 
性 求解 器 以 及 线性 求解 器 结构 . 文件 gcge_para.h 中 包含 程序 运行 过 程 中 所 需要 设置 的 
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各 种 参数 , 用 户 可 以 设置 其 中 的 参数 来 对 算法 运行 过 程 进行 控制 . 

第 三 层 的 头 文件 为 gcge_workspace.h, 该 文件 中 定义 了 GCGE_WORKSPACE 的 结 
W, 其 中 包含 了 算法 执行 过 程 中 所 需要 的 工作 空间 . 基于 gcge_ops.h 中 提供 的 创建 与 销 
毁 向 量 的 操作 , 以 及 gcge_para.h 中 的 参数 , 用 户 可 以 对 GCGE_WORKSPACE 中 的 工作 
空间 进行 创建 或 销毁 . 

依赖 第 二 、 三 层 的 头 文件 , 创建 了 第 四 层 gcge_rayleighritz.h, gege_orthogonalize.h, 
gcge_xpw.h 这 三 个 头 文件 , 其 中 包含 了 算法 运行 过 程 中 具体 用 到 的 函数 . 

依赖 于 前 四 层 头 文件 , 创建 了 第 五 层 的 文件 gcge_eigsol.h, 其 中 包含 具体 实现 GCG 算 
法 求解 特征 值 问题 的 主体 函数 . 

第 六 层 包 含 gcge_solver.h 文 件 , 其 中 定义 了 用 户 可 直接 使 用 的 GCGE_SOLVER 绪 
构 体 、 特 征 值 求解 函数 的 调用 接口 以 及 设置 用 户 可 调 参数 的 一 些 函 数 定义 . 


= 


2.1.2 SLEPc 格 式 的 调用 程序 结构 


SLEPc 格 式 的 调用 方式 如 图 1 所 示 . 其 中 右 侧 的 函数 以 及 中 间 部 分 的 函数 名 中 不 
包含 SLEPC 的 都 是 内 核 部 分 的 头 文件 gcge.h 中 给 出 的 , 中 间 部 分 函数 名 中 带 有 SLEPC 的 
是 app 部 分 的 头 文件 gcge_app_slepc.h 中 给 出 的 接口 函数 , 用 户 只 需 在 main 函 数 中 调用 
中 间 部 分 的 函数 就 可 以 进行 特征 值 问题 的 计算 . 

具体 调用 时 , 首先 读 入 矩阵 , 使 用 GCGE SOLVER SLEPC Create 创建 SLEPc 格 
式 的 求解 器 GCGE SOLVER, 设置 参数 , 然后 调用 GCGE SOLVER_SLEPC_Setup 函 
数组 装 之 后 , 就 可 以 调用 GCGE SOLVER Solve 函数 来 进行 特征 值 问题 的 求解 . 求解 结 
束 后 需要 调用 GCGE SOLVER SLEPO _ Free 函数 释放 求解 器 和 所 求 特征 对 所 占用 的 
内 存 空间 . 另外 , 用 户 可 以 通过 设置 打印 参数 来 确定 程序 运行 所 需要 输出 的 信息 , 如 各 部 
分 时 间 统 计 信息 与 特征 值 收 敛 情况 等 . 


SLEPC_ReadMatrixBinary 


GCGE_SOLVER_SLEPC_Create 


GCGE_PARA_SetFromCommandLine 
GCGE SOLVER SLEPC Setup 


GCGE SOLVER Solve 


| 


GCGE SOLVER Create 


GCGE SOLVER SetSLEPCOps 


GCGE SOLVER Setup 
GCGE EigenSolver 


GCGE SOLVER Free 


GCGE SOLVER SLEPC Free 


图 1: SLEPc-GCGE 调 用 方式 . 


3 算法 与 实现 过 程 的 优化 11 


3 ”算法 与 实现 过 程 的 优化 


本 节 基 于 对 GCG 算 法 性 质 的 分 析 与 理解 , 介绍 在 算法 实现 过 程 中 所 采用 的 一 些 用 来 
提高 算法 执行 效率 的 方法 , 主要 包括 高 效 实现 向 量 组 的 线性 组 合 、Rayleigh-Ritz 问 题 的 
求解 、 正 交 化 过 程 等 的 讨论 . 为 了 方便 描述 和 理解 , 这 里 采用 MATLAB 风 格 的 符号 来 进 
行 描述 . 


3.1 ”充分 使 用 BLAS-3 


在 算法 的 执行 过 程 中 , 我 们 充分 应 用 了 Level-3 的 BLAS 对 计算 过 程 进行 加 速 . 也 就 
是 对 于 串 行 计算 部 分 和 每 个 进程 内 部 的 计算 尽量 采用 和 矩阵 与 矩阵 相 乘 的 操作 . 对 于 
并 行 部 分 , 算法 中 用 到 的 所 有 分 布 式 存储 的 向 量 都 以 向 量 组 的 形式 出 现 . 只 要 用 户 的 
向 量 组 结构 支持 使 用 BLAS-3, 就 可 以 在 GCGE 中 直接 使 用 ， 比如 SLEPc 软 件 包 中 的 BV 
(BasisVectors) 结 构 , 其 局 部 问 量 数据 是 连续 存储 的 , 且 相 应 的 向 量 组 操作 使 用 了 BLAS- 
3 进行 加 速 , 这 个 BV 结构 就 可 以 直接 在 GCGE 软 件 包 中 使 用 . 目前 已 实现 的 SLEPc 接 口 
就 使 用 了 这 样 的 BV 结构 , 因此 在 已 实现 的 调用 接口 中 , 对 于 相同 矩阵 的 特征 值 进行 计算 
时 , SLEPc 接 口 的 计算 效率 是 最 高 的 . 

另外 对 于 用 户 只 提供 单 向 量 局 部 内 积 的 情况 , 我 们 的 多 向 量 内 积 操作 中 也 进行 了 统 
一 消息 传输 的 处 理 , 减少 了 并 行 时 进程 间 的 通讯 时 间 . 


3.2” 正 交 化 优化 


这 一 小 节 介 绍 软件 包 中 所 采用 的 正 交 化 方法 , 并 且 讨 论 如 何 提 高 向 量 组 正 交 化 的 计 
算 效率 与 并 行 效率 . 也 就 是 说 考虑 正 交 化 方法 的 时 候 需要 同时 考虑 正 交 化 过 程 的 效率 与 
稳定 性 , 尽量 在 这 两 个 方面 达到 一 个 好 的 平衡 , 以 使 最 后 的 计算 效率 达到 最 优 . 

除了 基本 的 Modified Gram-Schmidt [15] 正 交 化 方法 外 , 我 们 还 提供 了 Modified 块 
正 交 化 方法 ,Classical 块 正 交 化 方法 , 稳定 的 块 正 交 化 方法 , 以 及 多 重 正 交 化 方法 . 其 
中 Modified 块 正 交 化 方法 是 Modified Gram-Schmidt 正 交 化 方法 的 块 形式 , 两 者 可 以 达到 
相同 的 精度 , 但 是 块 正 交 化 的 方式 具有 更 高 的 计算 效率 与 并 行 效率 . Classical 块 正 交 化 方 
法 计算 量 小 , 并 行 效率 高 , 但 不 能 达到 必要 的 精度 而 使 得 特征 值 求解 算法 不 稳定 . 结合 以 
上 二 者 的 优势 得 到 了 稳定 的 块 正 交 化 方法 , 该 方法 可 以 得 到 较 高 的 精度 , 且 计算 效率 和 
并 行 效率 较 高 , 多 重 正 交 化 方法 几乎 可 以 与 Modified Gram-Schmidt 正 交 化 方法 达到 同 
样 的 精度 , 但 多 重 正 交 化 方法 的 计算 效率 和 并 行 效率 更 高 [18|. 

大 部 分 情况 下 , 我 们 推荐 使 用 多 重 正 交 化 方法 (默认 方法 ), 这 种 方法 在 第 4 节 的 测 
试 中 都 能 达到 需要 的 精度 ; 如 果 精 度 要 求 很 低 且 想 要 进一步 提高 效率 , 可 以 选择 使 
用 Classical 块 正 交 化 方法 进行 测试 ; 如 果 精 度 要 求 非常 高 , 使 用 默认 的 正 交 化 方法 无 法 收 
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Be, 可 以 使 用 Modified 块 正 交 化 方法 进行 测试 . 本 小 节 提 供 的 正 交 化 方法 在 软件 包 中 都 已 
实现 , 可 以 通过 设置 相应 的 正 交 化 参数 来 进行 选择 , 以 达到 用 户 的 稳定 性 和 精度 的 要 求 . 


3.2.1 Modified 块 正 交 化 方法 


Modified 块 正 交 化 方法 的 实现 过 程 描述 如 下 : 给 定 一 个 向 量 组 V, V 中 含有 的 向 量 
WVG, 1)... V( end), VH iistart — 1 个 向 量 已 经 是 在 矩阵 甩 意 义 下 的 单位 正 交 基 了 ， 
WV = V(41: start — 1), V = V(:, start : end)， 正 交 化 过 程 是 把 V 中 后 面 的 问 量 
组 VW 对 前 面 的 向 量 组 Vi 进行 正 交 , 然后 再 对 人 自己 进行 单位 正 交 化 . 

Modified 块 正 交 化 方法 具体 的 计算 过 程 描 述 见 算法 4, 其 中 reorth_count 表 示 重 正 交 
化 次 数 . 


Algorithm 4 Modified 块 正 交 化 方法 
1: EV, = B- Vi. 
2: for reorth count — 1: 2 do 
3: for i = 1 : start — 1 do > 去掉 他 中 人 方向 的 分 量 


4: 计算 prod = Va, (:, i)” V2. 

5: WHV = Vo — V(:, i)prod. 

6: end for 

7: for i = start : end do > Vo E £r QE SE 5 
8: 计算 prod = V(:, 1)? BV (i : end). 

9: 计算 V(:, 引 的 范 数 norm 一 Vprod(1). 

10: 如 果 norm > 7, 计算 V(: + 1: end) 2 V(:,7 +1: end) — 人 Vid) 

11: 否则 , 拷贝 V(:,i) =V(:,end — 1), 同时 令 end = end — 1, 4i — i — 1. 

12: end for 

13: end for 


正 交 化 算法 4 需要 的 临时 内 存 空 间 为 start 个 向 量 空间 与 end-start 个 浮 点 数 空 间 . 
观察 上 面 的 正 交 化 方法 可 以 发 现 里 面 的 for 循 环 步 对 每 个 向 量 进行 处 理 的 时 候 都 需要 进 
行 全 局 通讯 , 这 是 一 个 不 可 忽视 的 时 间 开 销 . 

Modified 块 正 交 化 方法 本 质 上 是 Modified Gram-Schmidt 正 交 化 方法 的 块 操作 . 
在 Modified Gram-Schmidt 正 交 化 方法 中 , 每 次 取 一 个 要 进行 正 交 化 的 向 量 , 减 去 在 
它 前 面 已 经 正 交 归 一 的 向 量 方向 上 的 分 量 . 为 了 不 损失 正 交 性 , 依次 减 去 每 个 向量 方向 
上 的 分 量 . Modifed 块 正 交 化 方法 则 是 同时 取 所 有 要 进行 正 交 化 的 向 量 , 依次 减 去 每 个 
已 经 正 交 归 一 向 量 方向 上 的 分 量 . Classical Gram-Schmidt 正 交 化 方法 相当 于 对 已 经 正 
交 归 一 的 向 量 组 做 一 个 块 处 理 , 这 样 可 以 使 计算 效率 变 好 , 但 往往 会 造成 数值 不 稳定 ; 
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Modified 块 正 交 化 方法 则 是 反 过 来 对 还 未 完成 正 交 归 一 的 向 量 组 做 块 处 理 , 这 样 既 可 以 
提高 求解 效率 , 又 能 保证 数值 稳定 性 . 


3.2.2 ”Classical 块 正 交 化 


在 分 布 式 并 行 计算 环境 下 , 使 用 Classical 块 正 交 化 方法 需要 进行 全 局 通讯 的 次 数 较 
少 , 计算 效率 高 . 但 由 于 数值 不 稳定 的 特点 (会 有 正 交 性 损失 ), 在 计算 精度 要 求 较 高 的 情 
况 下 往往 无 法 使 用 [15]. 

WV = V(:,1: start — 1), V = V(:, start : end). 算法 过 程 如 下 : 


Algorithm 5 Classical 块 正 交 化 方法 
: 计算 Ww = V - Vi - (VF BVÀ). > AHVE V A mkII RE 
计算 M = VI BV». > XV H ARTEZ 
: 求解 小 规模 特征 值 问题 MC = C9 的 全 部 特征 对 , 且 按 特征 值 从 小 到 大 排列 . 
: 初始 化 start V> = 1, 记 length_ 人 2 一 end 一 start +1. 
for i = 1 : length V5 do 

WRO; > 7, Ch: i) = ve; C. i). 

否则 , start Vo = start Vo- 1. 


: end for 
: EV = Vo: C(:, start_V2:length_V9). 


在 Classical 块 正 交 化 方法 中 , 去 掉 亿 中 仍 方向 分 量 的 方法 本 质 上 与 Classical Gram- 
Schmidt 正 交 化 方法 相同 , 会 有 一 定 的 正 交 性 损失 . 我 们 前 面 提 到 Classical Gram- 
Schmidt 正 交 化 方法 相当 于 对 已 经 正 交 归 一 的 向 量 组 做 一 个 块 处 理 ，Classical 块 正 交 
化 方法 是 既 对 已 经 正 交 归 一 的 向 量 组 Vi 做 块 处 理 , 也 对 未 完成 正 交 归 一 的 向 量 组 Ww 做 块 
处 理 , 因此 并 行 效率 可 以 进一步 得 到 提高 . 


3.2.3 ”稳定 的 块 正 交 化 方法 


为 了 既 提 高 正 交 化 的 计算 效率 , 又 能 保证 稳定 性 , 我 们 可 以 将 上 面 两 种 方式 进行 一 
定 的 结合 来 得 到 下 面 的 正 交 化 方法 . 

在 算法 6 的 第 二 步 中 , 对 VW 进行 正 交 归 一 化 的 时 候 采 用 的 是 与 Modified 块 正 交 化 一 
样 的 方式 而 使 得 Ww 具有 很 好 的 正 交 性 . 如 果 产 生 太 的 时 候 也 是 采用 Modified 块 正 交 化 的 
方式 而 导致 具有 很 好 的 正 交 性 , 这 样 也 会 使 得 在 进行 第 一 步 ( 减 去 WW 在 Wi 方向 上 的 分 
量 ) 计 算 的 时 候 具 有 很 好 的 数值 稳定 性 . 这 就 是 上 面 的 正 交 化 方法 虽然 在 第 一 步 中 采用 
了 块 处 理 的 形式 依然 具有 良好 的 稳定 性 的 原因 . 
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Algorithm 6 稳定 的 块 正 交 化 方法 

WHO = Vf BV. > EHV EV i RG 

: WHY, = V- VC. 

: 计算 矩阵 Frobenious 范 数 : Inner = ||C||p. 

: 如 果 Inner 大 于 正 交 的 靖 值 , 则 重复 上 面 三 个 步骤 的 计算 直到 Inner 小 于 阔 值 或 达到 
最 大 重复 次 数 . 

5: for i = start : end do > XIV; E AHIT EZ 

6: 计算 prod = V(:, 1)? BV (i : end). 

n 计算 V(:, 引 的 范 数 norm = Vprod(1). 

8 

9 


e w Ņ H 


如 果 norm > 7, 计算 V(: i +1: end) 2 V(5i-- 1: end) — prod(2iend) -V(:,2); 


否则 , 令 V(:,i) 2 V(: end), 同时 令 end = ena—1, i— i— 1. 
10: end for 


3.2.4 多重 正 交 化 


在 正 交 化 优化 过 程 中 可 以 发 现 尽量 对 已 收敛 的 向 量 组 进行 统一 正 交 化 ( 块 的 形式 ) 可 
以 减少 进程 间 消 息 发 送 的 次 数 , 同时 可 以 将 部 分 BLAS-2 的 操作 转化 为 BLAS-3 的 操作 . 
因此 我 们 考虑 使 用 二 分 递归 的 方式 来 尽量 对 已 收敛 的 向 量 进行 统一 的 正 交 化 处 理 . 这 
种 方法 可 以 将 几乎 全 部 的 BLAS-2 的 操作 转化 为 BLAS-3 的 操作 , 且 进 程 间 消 息 传输 的 次 
数 nc 与 向 量 组 中 的 向 量 个 数 nv 的 关系 为 : nc = O(nv). 由 于 这 里 设计 正 交 化 方法 的 思想 
来 源 于 多 重 网 格 和 迭代 算法 的 启发 , 所 以 我 们 将 这 种 方法 称 为 多 重 正 交 化 方法 . 后 来 我 们 
在 文献 综述 时 发 现 文章 [18] 也 讨论 了 类 似 的 算法 . 

考虑 对 问 量 组 V 中 的 第 start 位 置 到 end 位 置 的 向 量 进 行 正 交 化 , 算法 过 程 见 算法 7， 
其 中 size VV 表示 V 中 的 向 量 个 数 , [length/2| 表 示 不 小 于 length/2 的 最 小 的 整数 . 

算法 7 充分 利用 了 稳定 的 块 正 交 化 方法 的 优点 , 尽量 减少 要 进行 自身 正 交 化 向 量 组 
的 个 数 (只 对 单个 向 量 进行 正 交 归 一 化 ). 这 样 的 正 交 化 方式 除了 对 单个 向 量 进行 归 一 化 
外 其 它 的 部 分 都 可 以 用 BLAS-3 进 行 计算 , 从 而 可 以 充分 提高 计算 效率 和 并 行 效率 . 


3.3 ”减少 计算 子 空间 特征 值 的 个 数 与 向 量 组 线性 组 合 的 次 数 

本 小 节 将 目光 转向 如 何 高 效 地 计算 Rayleigh-Ritz 问 题 以 得 到 新 的 特征 向 量 逼 近 . 这 
个 过 程 包括 串 行 求解 一 个 小 规模 的 特征 值 问题 , 应 该 尽量 减少 其 在 整个 计算 时 间 中 所 占 
的 比重 . 求解 Rayleigh-Ritz 问 题 的 计算 过 程 包括 : 计算 小 规模 矩阵 4 = VT AV, 然后 求解 
标准 特征 值 问题 4C = CA, 最 后 将 计算 出 来 的 特征 向 量 C 返 回 到 长 向 量 组 V 进 行 向 量 的 
线性 组 合 以 得 到 新 的 特征 向 量 逼 近 X 与 向 量 组 已 . 
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Algorithm 7 多 重 正 交 化 Multi_Orth(V,， start, end) 
1: 计算 size_ V=end—start+1. 
2: Smid = start + [size_V/2] — 1, old_mid = mid. 


3: if size_V == 1 then 


4 计算 norm = ||V(:, start)||2. 

5 if norm > 7 then 

6: 计算 V(:, start) = V(:, start) /norm. 

7 else 

8 Send = start. 

9 end if 

10: else 

11: ， 调 用 Multi_ Orth(V,start,mid), 返回 V 和 mid. 


12: if mid < old_mid+1 then 


13: 记 n zero — old mid — mid, 更 新 end = end — n zero. 
14: 46 LV(:,mid: old_mid) = V(:, end : end-n zero). 
15: end if 


16: WY, = V(:, start : mid),Vo = V(:,mid +1: end). 

17: ”重复 计算 Ww = Vo — VVE - Va) BBA BITE SC EA BOR, 最 多 计算 3 次 . 
18: ”调用 Multi_ Orth(V,mid- 1, end). 

19: end if 

20: JKIBIV Mend. 


3.3.1 “计算 小 规模 矩阵 


在 求解 Rayleigh-Ritz 问 题 的 时 候 , 我 们 需要 先 把 相应 的 小 规模 矩阵 VTAV 和 VTBV 计 
算出 来 . 由 于 在 对 V 进 行 正 交 化 的 时 候 采 用 的 是 矩阵 B 的 内 积 , VTBV 就 是 一 个 单位 矩 
阵 , 所 以 不 需要 计算 . 那么 计算 小 规模 矩阵 的 过 程 主要 是 形成 矩阵 4 = VT AV. 

我 们 知道 V 的 形式 为 V = [X, P,W]. 为 了 与 前 面 算法 的 描述 相对 应 , 把 X 表 示 成 两 
部 分 X = [Xo, X], 其 中 Xo 表 示 已 收敛 的 特征 向 量 组 . 本 小 节 中 我 们 均 假 设 已 有 4 个 特征 
对 收敛 (0 < L < nev), Xi 表示 未 收敛 的 近似 特征 向 量 组 .由 算法 的 过 程 可 以 知道 , 这 里 
的 X 是 求解 上 一 次 迭代 的 Rayleigh-Ritz 问 题 得 到 的 , BUX = VC, 其 中 C 满 足 上 一 次 迭代 
中 的 Rayleigh-Ritz 问 题 : 


Agia? = CA, (1) 


其 中 C 是 一 个 正 交 和 矩阵 , AZ&— P XLBABER, Agate CAR Ee CD RE. 48 
阵 C 中 的 每 一 列 都 是 矩阵 4ua 的 一 个 特征 向 量 , 其 相应 的 特征 值 为 对 角 阵 A 相 应 位 置 的 对 
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角 元 素 . 由 此 可 以 知道 矩阵 4ua 有 如 下 的 谱 分 解 : 
Asta = CAC". (2) 


我 们 可 以 把 特征 向 量 组 C 分 解 成 两 部 分 C = [C,, CH, 其 中 0, 与 向 量 组 X 相 对应, 即 X = 
VC,. 与 X 的 结构 X = [Xo, Xi 相对 应 , 我 们 也 把 C, 分 解 成 Cs = [Co, Ci]. 
为 了 继续 进行 GCG 人 迭代 , 需要 形成 下 一 次 迭代 Rayleigh-Ritz 问 题 的 矩阵 4 如 下 : 
XTAX XTAP XTAW 
Anew =| P'AX PTAP PTAW |. (3) 
WTAX WTAP WTAW 
对 于 X74X 部 分 的 矩阵 块 有 如 下 的 性 质 


XTAX QV AMD SO Ag, = CICACTO, = 人。 (4) 


其 中 A。 表 示 在 对 角 阵 A 中 与 向 量 组 X 相 对 应 的 对 角 元 素 . 由 此 可 以 知道 我 们 并 不 需要 去 
直接 计算 矩阵 元 素 X7AX, MR REE RARE HIR AEA BA Anew P RIT. 

下 面 来 考虑 如 何 计算 矩阵 X74P 和 PT4P. 为 此 我 们 先 来 考虑 如 何 生 成 向 量 组 P. 
由 GCG 算 法 的 定义 可 以 知道 P = VC,, 其 中 0, 表示 与 向 量 组 P 相 对 应 的 系数 和 矩阵 0,. 为 
了 得 到 系数 矩阵 Cw 我 们 把 向 量 组 Ci 分 解 成 如 下 的 形式 : 


en eal 
Cis 
如 此 首先 可 以 设置 向 量 组 C, 如 下 : 
> 0 
Cp = : 
p Ga, 


接 下 来 对 小 规模 向 量 组 [Cu, C1, C, ETE L2 IE E46. 注意 这 里 的 [Ca, Ch] AA 8 IE IE HE, 
所 以 只 需要 把 C, 对 前 面 的 [Ca C1] 进 行 正 交 化 得 到 CG,, 然后 再 对 OG, 自 身 进行 正 交 化 即 可 . 
设 经 过 正 交 化 之 后 得 到 的 向 量 组 记 为 [Co, C1, Cp]. 注意 这 里 的 正 交 化 是 小 规模 向 量 组 的 
正 交 化 , 并 不 需要 对 长 向 量 组 [X, P| 直接 进行 正 交 化 . 

由 上 面 的 讨论 可 以 得 到 下 面 计 算 X+AP 和 PIAP 的 方式 : 


XAP = CTVTAVO, = CT Anal, = CTCACTC, = CTIC,, C2]A[C,, C41]7 C, 
= OF (CA0; + Cr Az (Or) Or = Cz Cs AC; Cp = 0, (5) 


和 


PT AP = C7 Ajay. (6) 
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最 后 由 (3), (4), (5) 和 (6) 可 以 得 到 矩阵 4ne。 的 结构 如 下 


Di 0 ay 
Anew = 0 a a2 , (7) 


T T 


其 中 oa = CT Aga4C,, a2 = WT AW, a, = XTAWAllag = PTAW. 也 就 是 说 这 里 的 aj， 
a 和 ao 是 需要 真正 针对 大 规模 向 量 组 进行 计算 的 , Di: 直接 就 是 对 角 和 矩阵 A 中 与 向 量 
组 X 相 对 应 的 部 分 A。, ai 可 以 通过 小 规模 矩阵 的 计算 C7 A aC pF BI. 

因为 Xo 己 经 达到 收敛 要 求 , 我 们 可 以 近似 地 认为 AXo = BXoAo, 于 是 近似 地 
&8 X1 AW = AoXLBW =0. 那么 (7) 近似 地 有 下 面 的 形式 


Ao 0 0 0 

= 0 A, O a 

Anew = i a ; (8) 
0 0 Q1 a2 


LT T 
0 aj a) O02 


FEHB Ag Be Cul BC ENT ARP IE TER | BOOGT £8 RB, A EAN A CS EE IE TECTA JS BE F8 AB E, 
ài = XL AW. 由 (8) 可 知 , ERE AL AE RIZR. 因此 只 需 对 矩阵 


的 特征 对 进行 求解 , 即 求 解 4swC = CA, 那么 未 sw 的 相应 的 特征 向 量 为 


I; 0 
(^2) " 


本 小 节 介 绍 的 方法 与 加 中 的 方式 基本 类 似 , 但 是 计算 aa 的 部 分 不 一 样 . 这 样 做 的 
原因 是 我 们 可 以 不 用 计算 和 矩阵 Awos 的 所 有 特征 对 , 而 只 需 计 算 Cs 部 分 即 可 . 且 使 用 优 
化 (8) 时 可 以 使 实际 要 求解 特征 值 问 题 的 矩阵 规模 更 小 , 这 样 就 为 下 一 步 减少 计算 小 规模 
FEE AE IRISUA, C = CA 的 时 间 提 供 了 条 件 . 


3.3.2 ”使 用 dsyevx 代 兰 dsyev 来 求解 部 分 特征 值 


在 上 一 小 节 计 算出 Rayleigh-Ritz 特 征 值 问 题 的 矩阵 4 之 后 , 接 下 来 就 是 要 求解 特 
征 值 问 题 4j6sC = CA. 最 常用 的 方式 就 是 调用 LAPACK |1] 中 的 函数 dsyev 来 求解 所 有 特 
征 对 . 为 了 尽量 减少 计算 量 , 我 们 只 求解 子 空 间 矩 阵 Apes 的 第 L 二 1 到 第 nev 个 特征 值 和 相 
应 的 特征 向 量 C,, 具体 是 调用 dsyevx 来 进行 计算 . 当然 在 下 一 小 节 的 讨论 中 , 我 们 会 把 这 
一 部 分 的 计算 量 进一步 地 减少 . 
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3.3.3 ”进一步 减少 计算 子 空间 特征 值 的 个 数 与 向 量 组 线性 组 合 的 次 数 
在 计算 Rayleigh-Ritz 问 题 的 过 程 中 , 我 们 需要 用 串 行 的 方式 求解 如 下 特征 值 问题 : 


Aust ON (11) 


VE FH t] de SERE HH, 不 同 的 特征 对 是 依次 收敛 的 . 由 于 V 中 会 保存 已 收敛 的 特征 向 量 ， 
所 以 对 于 特征 值 问题 (11), 相应 的 单位 向 量 (0,.… 0, 1, 0, .…0) 就 是 其 相应 于 已 收敛 特征 对 
的 一 个 特征 向 量 , 其 中 1 所 在 的 位 置 即 为 这 个 已 收敛 特征 对 的 编号 . 也 就 是 说 在 求解 特征 
值 问题 (11) 的 时 候 , 对 已 经 收敛 的 特征 向 量 不 需要 再 进行 求解 了 , 而 只 需要 求解 还 未 收敛 
的 那些 特征 向 量 , 即 与 X1 相 对 应 的 特征 对 . 这 样 就 减少 了 调用 dsyevx 计 算 特征 对 的 个 数 
而 进一步 减少 计算 时 间 . 

更 进一步 , 在 已 经 有 特征 向 量 收敛 的 情况 下 , 当 进 行 小 规模 向 量 组 正 交 化 ( 计 
算 [C。, Cy]) 的 时 候 可 以 进一步 减少 计算 量 . 设 已 有 4 个 特征 对 收敛 , 即 特 征 值 问题 (11) 的 
前 6 个 特征 向 量 为 [1,0]7, 所 以 只 需要 将 新 得 到 的 子 空间 向 量 组 [C,, Co 每 列 的 前 个 元 素 
赋 为 0, 就 可 使 其 与 前 (个 特征 向 量 正 交 . 然后 再 对 这 些 子 空间 特征 向 量 组 的 非 零 部 分 进 
行 正 交 化 就 可 以 得 到 正 交 向 量 组 [C,, C]. 

为 了 方便 理解, 我 们 通过 下 面 的 具体 过 程 来 阐述 上 面 的 方法 .首先 通过 算法 的 构造 
过 程 可 以 知道 矩阵 未 。 具 有 如 下 的 结构 ; 


Di 0 ay 
Anew = 0 [941 Q2 


T T 


在 迭代 过 程 中 , 对 于 第 i 个 已 收敛 的 特征 对 , 其 特征 向 量 已 达到 收敛 准则 , 可 以 不 必 更 
新 , 因此 可 将 其 子 空间 特征 同 量 记 为 e;， 设 已 有 l 个 连续 的 特征 对 收敛 , 那 只 需要 再 求 
fH Anew HC + 1 到 dimX 个 特征 对 . 这 里 dimxX 表 示 向 量 组 [Xo, Xi] 中 的 向 量 个 数 . 那么 可 以 
知道 4w6s 的 前 dimx 个 特征 向 量 如 下 


要 使 XX 的 各 列 均 正 交 , 首先 将 Cis 设 为 0, 再 对 Cs 进行 正 交 化 以 得 到 新 的 向 量 组 C2. 记 正 
交 化 后 的 特征 向 量 为 
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因为 < dimx, 所 以 Cs 可 以 写成 下 面 的 形式 

lL 0 
0 
O° Gy 
0 qm 


KB Cu, Cop 和 Cou 分 别 与 向 量 组 Y 中 的 X, P 和 W 部 分 相对 应 . 这 样 就 可 以 构造 出 如 下 的 


与 长 向 量 组 [X, PI 所 对 应 的 系数 向 量 组 


= 0 
Cz, Cpl] E 0 
0 Co 

4 需要 对 下 面 的 子 块 进行 正 交 化 即 可 


由 此 可 以 知道 要 对 P 向 量 组 部 分 进行 正 交 化 , 只 


对 上 面 的 向 量 组 进行 正 交 化 之 后 得 到 如 下 的 向 量 组 


然后 利用 存储 在 V 中 的 基 向 量 组 进行 线性 组 合 以 得 到 新 的 向 量 组 [X, P] = VI 
据 上 面 给 出 的 系数 矩 阵 [G,, C], 对 于 已 经 收敛 的 特征 方向 X(:, 1 : 0 不 需要 进 


且 X(:,《 十 1 :dimX) 的 更 新 采用 如 下 方式 


Cs 
X(:,€+1:dimx) =V(:,€+1:dimXPW)-| Cæ |, 


问 量 组 P 的 更 新 方式 如 下 


P=V(:,€+1: dimXPW) - ma Ae 
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其 中 qimXPW 表 示 基 向 量 组 [X, PW] 中 的 向 量 个 数 . 

如 果 一 个 已 收敛 的 特征 值 是 重 的 , 且 它 的 特征 空间 还 未 完全 收敛 , 那么 需要 对 重 特 
征 值 进行 一 定 的 处 理 . 在 检查 未 收敛 特征 对 的 收敛 性 时 , 我 们 同时 对 相应 特征 值 的 重 
数 进行 检查 .如 果 这 个 特征 值 与 前 一 个 特征 值 的 相对 误差 小 于 一 个 阔 值 (这 里 默认 的 阐 
值 1e-6), 那 就 认为 这 个 特征 值 与 前 一 个 特征 值 是 重 的 . 如 果 当 前 检查 的 特征 对 已 收敛 ， 
那么 就 检查 其 与 前 一 个 特征 值 是 否 是 重 特征 值 , 如 果 是 重 的 且 前 一 个 特征 对 被 标记 为 未 
收敛 , 那么 就 将 当前 的 特征 对 也 标记 为 未 收敛 ; 如 果 当 前 检查 的 特征 对 未 收敛 , 也 检查 其 
与 前 一 个 特征 值 是 否 是 重 特征 值 , 如 果 是 重 的 且 前 一 个 特征 对 被 标记 为 已 收敛 , 那么 就 
将 前 一 个 特征 对 也 标记 为 未 收敛 . 检查 收敛 性 直到 连续 未 收敛 的 特征 对 个 数 达 到 一 定 的 
上 限 (默认 为 10), 之 后 的 特征 对 都 标记 为 未 收敛 , 或 者 到 最 后 一 个 特征 对 为 止 . 记 从 第 1 个 
特征 对 开始 连续 收敛 的 最 后 一 个 特征 对 的 编号 为 {is, 则 求解 子 空间 矩阵 特征 值 问 题 (11) 
时 要 从 第 pain 十 1 个 开始 计算 , PSW 向 量 组 的 构造 是 相应 于 XX 中 标记 为 未 收 化 的 向 量 进 
行 . 


3.4 ”计算 子 空间 特征 值 问 题 的 优化 
3.4.1 ”只 让 0 号 进程 计算 子 空间 特征 值 问题 

由 于 各 进程 计算 的 是 相同 子 空间 上 的 Rayleigh-Ritz 特 征 值 问题 , 且 计 算 完毕 后 需要 
由 0 号 进程 把 相应 的 特征 对 广播 到 各 个 进程 以 保持 所 有 进程 的 数据 统一 , 在 广播 前 可 能 
会 有 进程 等 待 的 情况 , 所 以 只 由 0 号 进程 进行 计算 , 可 以 在 一 定 程度 上 减少 这 部 分 的 计算 
时 间 . 


3.4.2 ”每 个 进程 计算 一 部 分 特征 值 , 再 统一 到 所 有 进程 


当 求 解 特征 值 个 数 较 多 时 , 由 于 计算 子 空间 特征 值 问题 的 时 间 与 求解 特征 值 的 个 数 
是 超 线性 关系 , 这 部 分 的 计算 时 间 占 比 会 变 得 很 高 . 为 了 减少 计算 时 间 , 我 们 将 特征 值 问 
题 (11) 平 均 分 配 到 所 有 进程 进行 求解 , 每 个 进程 只 计算 少量 的 几 个 特征 对 , 这 样 就 可 以 大 
大 地 减少 求解 子 空间 特征 值 问 题 (11) 的 计算 时 间 . 计算 结束 后 , 使 用 MPI_Allgatherv 将 
所 有 进程 的 特征 对 信息 进行 统一 . 在 进程 数 不 是 那么 多 的 时 候 , 计算 特征 对 的 时 间 相 
比 MPI_Allgatherv 的 时 间 要 多 很 多 , 此 时 采用 这 样 的 方式 就 可 以 提高 计算 效率 . 为 了 使 
所 有 进程 的 子 空间 特征 对 信息 统一 , 我 们 需要 对 它们 进行 广播 . 这 部 分 的 消息 传输 是 必 
不 可 少 的 , 当 进 程 数 很 多 的 时 候 , 可 以 考虑 只 让 其 中 一 部 分 的 进程 进行 计算 从 而 减少 消 
息 传 递 所 占 的 时 间 比 重 . 为 了 避免 一 些 进程 中 计算 的 特征 值 个 数 太 少 又 需要 其 对 所 有 进 
程 发 送 它 所 计算 的 特征 值 , 我 们 设置 让 每 个 进程 至 少 计 算 一 定 个 数 的 特征 值 , 默认 每 个 
进程 至 少 计算 10 个 . 
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3.5 ”使 用 MKL 进 行 优化 


为 了 进一步 减少 计算 时 间 , 当 使 用 Inte] 编 译 器 时 可 以 调用 mkl 的 BLAS 与 LAPACK 库 . 
一 般 地 , 如 果 计 算 设 备 上 用 的 是 Intel 的 CPU 并 且 同 时 有 mkl 的 BLAS 和 LAPACK 库 的 时 
候 , 建议 使 用 它们 . 


4 测试 结果 


以 下 进行 的 所 有 测试 的 计算 环境 为 中 国 科 学 院 数 学 与 系统 科学 研究 院 科 学 与 工程 
计算 国家 重点 实验 室 LSSC4 计算 集群 , 提交 任务 时 采用 的 是 big 队 列 . 每 个 节点 包括 2 颗 
主 频 为 2.3GHz 的 Intel Xeon Gold 6140 18 核 Purley 处 理 器 和 192GB 内 存 , 更 详细 的 介绍 可 
以 参看 http://lsec.cc.ac.cn/chinese/lsec/LSSC-IVintroduction.pdf. 

我 们 对 以 下 四 组 对 称 和 矩阵 进行 了 测试 . 其 中 前 3 个 矩阵 取 自 Suite Sparse Matrix 


表 1: 测试 矩阵 . 

编号 JERE 维 数 ” 非 零 元 个 数 
1 Andrews 60000 760154 
2 Ga3As3H12 61349 5970947 
3 Gal0As10H30 113081 6115633 

4 RITER 2736987 189967149 


Collectionl，Andrews 和 矩阵 的 第 一 个 特征 值 接近 于 0, Ga3As3H125 GalOASs10H305B BEY 
有 大 量 负 特征 值 , 这 三 个 矩阵 的 特征 值 分布 均 非常 稠密 ; 模 态 分 析 和 矩阵 是 从 某 建筑 体 的 
模 态 分 析 中 导出 的 广义 代数 特征 值 问题 矩阵 [6, 7, 19]. 

在 软件 包 GCGE 所 提供 的 接口 中 , 进行 了 向 量 组 计算 优化 的 SLEPc 接 口 是 计算 效率 
最 高 的 , 因此 以 下 的 数值 计算 测试 均 使 用 SLEPc 的 调用 接口 执行 . 为 了 进行 对 比 , 我 们 还 
使 用 SLEPc 中 的 LOBPCG 以 及 Jacobi-Davidson 解 法 器 进行 了 测试 . 


4.1 ”优化 策略 对 计算 效率 提升 的 测试 

在 这 一 人 小节, 我 们 用 实际 的 数值 算 例 来 测试 本 文中 介绍 的 算法 优化 的 效果 . 这 里 求 
解 的 特征 值 问题 中 的 和 矩阵 为 Andrews 和 矩阵 , 收敛 准则 为 1e-4, 在 单个 节点 上 用 24 个 进程 求 
解 前 600 个 最 小 端的 特征 值 和 相应 的 特征 向 量 . 表 2 展 示 了 使 用 文中 优化 过 程 所 带 来 时 间 
减少 的 情况 , 其 中 RR 过 程 指 Rayleigh-Ritz 过 程 的 计算 时 间 . 


详 见 https://sparse.tamu.edu 


mn 
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分 批 计 算 (120) 表 示 进 行 分 批 计 算 时 设置 参数 block_size 为 120, 根据 经 验 , 我 们 
设置 分 批 计算 时 参数 block_size 的 默认 值 为 nev/5. 这 种 方式 会 减少 每 次 迭代 中 向 量 
组 P 与 向 量 组 W 的 计算 量 , 同时 也 可 以 减 小 Rayleigh-Ritz 过 程 的 计算 规模 , 从 计算 结果 
可 以 看 到 , 计算 P, 计算 W 以 及 Rayleigh-Ritz 过 程 的 计算 时 间 均 有 减少 . 

向 量 运 算 优 化 表示 对 算法 中 同 量 的 操作 尽量 使 用 向 量 组 统一 计算 的 方式 , 这 样 可 以 
尽 可 能 地 使 用 BLAS-3 对 算法 进行 加 速 . 从 计算 结果 可 以 看 到 , 算法 各 部 分 的 计算 时 间 均 
有 明显 减少 , 尤其 是 计算 XX 与 P 的 部 分 , 这 是 因为 这 两 部 分 的 计算 量 主 要 是 进行 可 以 完全 
转化 为 BLAS-3 操 作 的 向 量 组 线性 组 合 的 操作 . 

在 正 交 化 优化 中 , W 的 正 交 化 采用 的 是 不 稳定 但 最 高 效 的 Classical 块 正 交 化 方法 , BI 
算法 5, P 的 正 交 化 采用 的 是 稳定 且 高 效 的 多 重 正 交 化 方法 , 即 算法 7. 因此 计算 P 和 W 两 
部 分 的 计算 时 间 有 明显 地 减少 . 

最 后 的 稠密 特征 值 优 化 , 是 在 Rayleigh-Ritz 过 程 中 使 用 不 同 进程 分 别 求解 一 部 分 
特征 对 , 再 统一 进行 MPI 的 全 收集 操作 . 这 使 得 Rayleigh-Ritz 过 程 的 计算 时 间 有 明显 地 
减少 . 


K 2: Andrews 和 矩阵 , 收敛 准则 取 1le-4, 使 用 不 同 优化 方法 的 时 间 对 比 ( 秒 ). 
优化 方式 计算 X 计算 已 计算 W RR 过 程 总 时 间 

原始 GCG 25.06 32.98 101.50 33.42 193.49 

分 批 计算 (120) 23.20 849 52.97 15.43 101.77 

向 量 运算 优化 251 2.05 50.81 10.19 65.98 

正 交 化 优化 1.42 0.84 16.74 9.79 29.20 

Bel a xp RE ETE DUI 143 0.83 16.94 413 28.82 


4.2 求解 特征 值 个 数 与 计算 时 间 关 系 的 测试 


为 了 测试 GCG 算 法 求解 特征 值 问题 的 时 间 与 所 求 特征 值 个 数 的 关系 , 我 们 使 用 一 
个 节点 上 的 24 个 进程 对 表 1 中 的 前 3 个 矩阵 进行 了 测试 . 分 别 计算 了 75-1200 个 特征 值 , 测 
试 中 均 使 用 了 绝对 残 差 (4z 一 和 Bzjlz/|z|l2) 作 为 收敛 准则 , 分 别 对 比 了 收敛 准则 取 le- 
4 和 1e-12 的 结果 . 这 里 测试 使 用 GCG 算 法 时 , 分 批 计算 的 参数 均 取 block_size = nev/5, 
且 求 解 线性 方程 组 时 均 使 用 10 次 CG 迭 代 , 使 用 的 正 交 化 方法 均 为 多 重 正 交 化 方法 (这 
种 正 交 化 方法 在 精度 、 稳 定性 和 可 扩展 性 方面 达到 了 很 好 的 平衡 ). 由 于 LOBPCG 方 法 
所 展现 出 来 的 良好 的 计算 效率 与 可 扩展 性 , 本 小 节 与 4.3 小 节 的 计算 结果 均 与 SLEPc 中 
的 LOBPCG 进 行 了 比较 . 
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4.2. Andrews 和 矩阵 


对 Andrews 和 矩阵 , 收敛 准则 取 1le-4 时 . 经 过 实际 的 数值 测试 , 我 们 发 现 LOBPCG 的 参 
Zi—eps lobpcg _ blocksize 设 置 为 nev/6, —eps lobpcg restartix E7J0.l1HJH] fixit 
算 效率 最 高 . 我 们 采用 这 样 的 参数 设置 所 得 到 的 计算 结果 来 进行 比较 . 且 在 此 之 后 的 参 
数 一 eps_lobpcg_restart 均 设置 为 0.1. 两 种 方法 的 计算 时 间 对 比 情况 见 表 3 与 图 2， 从 
计算 结果 的 对 比 可 以 发 现 GCGE 相 比 LOBPCG, 计算 效率 有 约 3-4 倍 的 提升 . 


K 3: Andrews 和 矩阵 , 收敛 准则 取 1e-4 时 两 种 方法 的 计算 时 间 对 比 ( 秒 ). 
nev GOGE LOBPCG 


75 1.98 10.12 
150 4.51 24.42 
300 . 10.59 60.95 


600 28.20 146.62 
1200 103.56 405.38 


CPU time for Andrews, tol = 1e—4, np=24 


erai slope=1.5 
=i GCGE time 
= LOBPCG(SLEPc) time 


CPU time (in seconds) 


Number of eigenvalues 


图 2: Andrews Ee, 收敛 准则 取 1le-4 时 两 种 方法 的 计算 时 间 对 比 . 


收敛 准则 取 1e-12 时 , LOBPCG 如 采 仍 使 用 收敛 准则 为 le-4 时 的 最 优 参 数 , 在 特 
征 值 个 数 分 别 为 150,600,1200 时 计算 失败 , 并 且 特 征 值 个 数 分 别 为 75,300 时 , 将 参 
Zi—eps lobpcg _blocksize 设 置 为 nev 或 nev/6 的 时 候 计 算 时 间 几 平一 致 . 因此 将 参 
数 一 eps_lobpcg_blocksize 设 置 为 nev 来 对 LOBPCG 进 行 测试 . 两 种 方法 计算 时 间 的 对 
比 情况 见 表 4 与 图 3. 由 其 中 的 数值 结果 可 以 发 现 , GCGE 计 算 更 稳定 , 且 计 算 效 率 相 对 
于 LOBPCG 提 升 约 6 倍 . 


4 测试 结果 24 


K 4: Andrews 和 矩阵 , 收敛 准则 取 1e-12 时 两 种 方法 的 计算 时 间 对 比 ( 秒 ). 
nev GCGE LOBPCG 
75 5.46 24.20 
150 11.54 62.75 
300 . 28.33 168.08 
600 70.63 512.89 
1200 249.04 1931.54 


CPU time for Andrews, tol = 1e—12, np-24 


LE slope=1.5 
10°} | =i GCGE time 
= LOBPCG(SLEPc) time 


CPU time (in seconds) 


Number of eigenvalues 


图 3: Andrews#E RE, 收敛 准则 取 1le-12 时 两 种 方法 的 计算 时 间 对 比 . 


4.2.2 Ga3As3H12 和 矩阵 

本 小 节 对 Ga3As3H12 和 矩阵 进行 测试 . 24 uic S VE MUI gr = le-4 时 ，LOBPCG 的 参 
数 -eps_1lobpcg_blocksize 设 置 为 nev/6, 相应 的 计算 时 间 对 比 情况 见 表 5 与 图 4， 由 
其 中 的 结果 我 们 可 以 发 现 GCGE 的 计算 效率 比 LOBPCG 有 约 2-3 倍 的 提升 . 


K 5: Ga3As3H12 和 矩阵 , 收敛 准则 取 1e-4 时 两 种 方法 的 计算 时 间 对 比 ( 秒 ). 
nev GCGE LOBPCG 


75 8.36 20.73 
150 14.91 41.88 
300 27.80 84.33 


600 57.52 160.32 
1200 148.84 318.99 


SUCHE 7r = le-12 时 , 使 用 LOBPCG 进 行 计算 的 时 候 往往 失败 或 者 在 一 定时 间 
内 无 法 收敛 , 计算 结果 对 比 情 况 见 表 6, 其 中 LOBPCG 的 参数 -eps_lobpcg_blocksize 
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, CPU time for Ga3As3H12, tol = 1e-4, np=24 
10 


Uer wet Slope-1.5 
—f-— GCGE time 
= LOBPCG(SLEPc) time 


CPU time (in seconds) 


Number of eigenvalues 


图 4: Ga3As3H12 和 矩阵 , 收敛 准则 取 1le-4 时 两 种 方法 的 计算 时 间 对 比 . 


设置 为 nev. 此 时 可 以 发 现 GCGE 有 明显 的 稳定 性 和 效率 的 优势 . 


K 6: Ga3As3H125B ME, 收敛 准则 取 1e-12 时 两 种 方法 的 计算 时 间 对 比 ( 秒 ). 


nev GCGE LOBPCG 
75 23.45 3000025 cat T 197 
150 40.78 Error in Lapack DSYGVD 127 
300 75.67 运行 超过 24 小 时 候 未 收敛 
600 149.93 运行 超过 24 小 时 候 未 收敛 
1200 368.65 运行 超过 24 小 时 候 未 收敛 


4.2.3 Gal10As10H305B[£ 

本 小 节 对 Gal0Asl10H30 和 矩阵 进行 测试 , 24 Ac ot HE Wr = 1e-4 时 ,， LOBPCG 的 参 
数 -eps_lobpcg_blocksize 设 置 为 nev， 两 种 方法 的 计算 时 间 对 比 情况 见 表 7 与 图 5. 
由 其 中 的 数值 结果 可 以 发 现 GCGE 的 计算 效率 相对 于 LOBPCG 有 约 2-3 倍 的 提高 . 


表 7: Gal0As10H30 和 矩阵 , 收敛 准则 取 1le-4 时 两 种 方法 的 计算 时 间 对 比 ( 秒 ). 
nev GCGE LOBPCG 
75 9.37 24.18 
150 20.07 52.36 
300 40.32 109.53 
600 78.72 231.73 
1200 205.24 559.18 
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CPU time for Ga10As10H30, tol = 1e-4, np-24 
10 r , 


eae slope=1.5 
== GCGE time 
—Ei- LOBPCG(SLEPo) time 


CPU time (in seconds) 


10° 10° 
Number of eigenvalues 


图 5: Gal10As10H30 和 矩阵 , 收敛 准则 取 1e-4 时 两 种 方法 的 计算 时 间 对 比 . 


当 收 敛 准则 为 + = 1e-12 时 , 使 用 LOBPCG 进 行 求解 往往 计算 失败 或 者 在 一 定时 间 内 
无 法 收敛 , 计算 结果 的 对 比 情况 见 表 8, 其 中 LOBPCG 的 结果 是 参数 -eps lobpcg blocksize 
设置 为 nev 的 时 候 得 到 的 . 同样 可 以 发 现 GCGE 有 明显 的 稳定 性 和 效率 上 的 优势 . 


表 8: Gal0As10H30 和 矩阵 , 收敛 准则 取 1e-12 时 两 种 方法 的 计算 时 间 对 比 ( 秒 ). 
nev GCGE LOBPCG 

75 25.17 Error in Lapack DSYGVD 78 
150 55.16 Error in Lapack DSYGVD 134 


300 110.07 运行 超过 24 小 时 未 收敛 
600 212.47 运行 超过 24 小 时 未 收敛 
1200 526.26 运行 超过 24 小 时 未 收敛 


4.3 ” 强 可 拓展 性 测试 

为 了 进行 算法 强 可 拓展 性 的 测试 , 我 们 对 Andrews 和 矩阵 ， Gal0As10H3048 84, 分 别 使 
用 1-32 个 进程 进行 测试 . 这 里 同样 将 GCGE 与 SLEPc 中 的 LOBPCG 进 行 对 比 , 具体 的 计 
算 结 果 见 图 6 与 图 7. 对 比 图 6 和 7 可 以 发 现 , GCGE 有 与 LOBPCG 几 平一 样 的 强 可 拓展 性 ， 
并 且 GCGE 具 有 更 好 的 计算 效率 . 


44 ”大 规模 计算 时 算法 的 稳定 性 与 计算 效率 测试 
这 一 小 节 考 虑 求解 茶 种 复杂 建筑 体 横 态 分 析 所 导致 的 特征 值 问题 . 模 态 分 析 问 题 的 
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CPU time for Andrews, tol = 1e-4, nev=600 


renee slope--0.77 
—i— GCGE time E 
= LOBPCG(SLEPc) time 


CPU time (in seconds) 


10° 10" 
Number of processers 


图 6: Andrews 和 矩阵 , 收敛 准则 取 1le-4, 与 LOBPCG 强 可 拓展 性 对 比 . 


GPU time for Ga10As10H30, tol = 1e-4, nev=600 
10 


bist Slope=-0.81 
一 此 一 GCGE time 
= LOBPCG(SLEPc) time 


CPU time (in seconds) 


10° 10" 
Number of processers 


图 7: Gal0As10H305R pE, 收敛 准则 取 1le-4, 与 LOBPCG 强 可 拓展 性 对 比 . 


力学 模型 如 下 
-GAu — (u+ G)V(V -u) = Au, inQ, (12) 
u = 0, onðQ, 
其 中 
» ae aaa hii 
= — o 。 一 — a, = 
= A jS dy 2 PITE Uy 
A 2 uz 


REGIA EB WR E (RERE), /为 材料 体积 模 量 , RE BLEUS (Lamé 
constant), 其 表达 式 如 下 : 
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其 中 五 为 弹性 体 的 杨 氏 模 量 , z 为 Poisson 比 . 由 于 弹性 体 可 能 由 多 种 材料 构成 , 所 以 拉 梅 
常数 在 区 域 的 不 同位 置 可 能 取 不 同 的 值 . 本 小 节 考 虑 的 力学 模型 有 如 下 的 应 用 特征 


表 9: 模 态 分 析 问 题 的 部 分 应 用 特征 . 
模型 Poisson 比 杨 氏 模 量 
某 建筑 体 模 态 分 析 问 题 ”0.15, 0.3 7.0e7, 2.48e7, 2.1e8 


本 小 节 所 测试 的 矩阵 是 使 用 表 9 中 的 参数 , 并 用 有 限 元 并 行 计算 框架 Panda [7, 19] 进 
行 离散 得 到 的 . 这 组 矩阵 是 对 称 正定 的 , 但 由 于 使 用 了 多 种 材料 , 且 网 格 剖 分 的 精细 程度 
不 同 , 导致 该 组 矩阵 具有 很 强 的 多 尺度 特性 . 这 就 使 得 矩阵 条 件数 很 大 , 对 解法 器 的 稳定 
性 要 求 很 高 , 给 特征 值 的 计算 带 来 了 很 大 的 挑战 . 表 10 中 给 出 了 这 组 矩阵 的 部 分 数学 性 
质 , RPA, B 分 别 表示 该 广义 特征 值 问 题 的 刚度 矩阵 与 质量 矩阵 .由 表 10 可 以 看 出 这 组 
和 矩阵 具有 很 强 的 多 尺度 特性 且 条 件数 很 大 . 


表 10: 模 态 分 析 问 题 所 产生 和 矩阵 的 性 质 . 


性 质 A B 

行 数 273,6987 273,6987 

非 零 元 个 数 1,8996,7149  1,8996,7149 

稀疏 度 69.41 69.41 

对 角 线 元 素 最 小 值 9.90e--08 3.88e-02 

对 角 线 元 素 最 大 值 3.45e+11 4.29e--01 

非 对 角 线 元 素 最 小 值 -1.67e+11 2.34e-03 
非 对 角 线 元 素 最 大 值 1.69e+11 1.60e--01 
非 对 角 线 元 素 最 小 绝对 值 2.91e-11 2.34e-03 
非 对 角 线 元 素 最 大 绝对 值 1.69e--11 1.60e--01 
最 小 对 角 占 有 率 1.32e-02 2.60e-02 

对 称 性 sym (1.0e+00) sym (1.0e-10) 

条 件数 8.57e--07 7.74e 十 03 


我 们 使 用 SLEPc 软 件 包 中 的 几 种 常用 的 算法 对 本 小 节 的 问题 进行 了 测试 . 使 
用 Krylov-Schur 方 法 求解 时 , 可 以 达到 要 求 的 计算 精度 , 且 收敛 速度 快 . 缺点 是 该 方法 需 
要 精确 求解 迭代 过 程 中 的 线性 方程 组 , 由 于 本 小 节 的 矩阵 的 条 件数 较 大 , 只 能 使 用 LU 分 
解 才 能 求解 , 因此 并 行 可 拓展 性 较 差 . LOBPCG 方 法 不 需要 精确 求解 线性 方程 组 , 可 以 
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有 较 好 的 并 行 可 拓展 性 , 但 该 方法 稳定 性 较 差 , 求解 本 小 节 的 矩阵 特征 值 问 题 时 总 是 求 
解 失败 . Jacobi-Davidson 方 法 不 需要 精确 求解 线性 方程 组 , 并 行 可 扩展 性 较 好 , 但 该 方法 
较为 复杂 , 算法 参数 很 多 , 难以 为 本 小 节 的 矩阵 找到 最 优 的 求解 参数 以 及 构造 良好 的 预 
AK TAB IE. 

本 小 节 使 用 2736987 阶 的 矩阵 [6, 了] 进行 测试 , 收敛 准则 为 相对 残 差 ||Ax 一 和 Bz||2/ (|A|: 
|zll2) < 5e-2, 使 用 72 到 576 个 进程 , 分 别 求解 最 小 端的 100 和 1000 个 特征 值 . 根据 上 面 的 
分 析 , 这 里 我 们 将 GCGE 与 SLEPc 中 的 Jacobi-Davidson 解 法 器 进行 对 比 . 使 用 GCGE 求 
解 时 , 用 500 次 CG 和 迭代 计算 全 . 经 过 多 次 试验 , 也 为 Jacobi-Davidson 方 法 选取 了 已 知 最 好 
的 参数 . 相应 的 对 比 结果 见 图 8. 从 图 中 可 以 发 现 GCGE 与 Jacobi-Davidson 均 有 较 好 的 强 
可 拓展 性 , 且 GCGE 的 计算 速度 相 比 Jacobi-Davidson 方 法 有 约 3-5 倍 的 提高 . 


CPU time for guangji, tol = 5e-2 
JD time, nev=1000 


GCGE time, nev=1000 
=EF UD time, nev=100 
-= GCGE time, nev-100 


ES 
O 
in 


CPU time (in seconds) 


10° 
Number of processers 


图 8: 2736987 HEAS 2) Hr FE BE, 收敛 准则 取 5e-2 时 的 计算 时 间 对 比 . 
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本 文 基于 阻尼 块 反 野 法 与 子 空间 投影 方法 的 组 合 构造 了 求解 特征 值 问题 的 广义 共 
He RE BGA, 并 构造 了 相应 的 算法 软件 包 GCGE. 该 软件 包 具 有 稳定 、 高 效 、 高 可 拓展 
性 的 特点 . 针对 不 同 的 问题 , 计算 效率 相 比 于 SLEPc 软 件 包 中 的 LOBPCG 以 及 Jacobi- 
Davidson 解 法 器 有 2-6 倍 的 效率 提升 . 且 该 软件 包 是 Matrix-Free 和 Vector-Free 的 , 因此 可 
以 有 更 广泛 的 应 用 , 未 来 可 以 依据 阻尼 块 反 窜 法 设计 求解 非 对 称 和 矩阵 特征 值 问 题 的 相应 
算法 . 
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